
C
C funcao para calculo de problemas de potencial bi-dimensionais pelo
C metodo dos elementos de contorno com elementos constantes.
C

	FUNCTION RESID(PAR,fvec,dudQOT,m,numpar,ldfjac,iflag,metodo)
	
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)


	INCLUDE 'param_dim.inc'

	
	CHARACTER*11 FILE_B
	CHARACTER*15 FILE_DAT,FILE_DCC,FILE_DEX,FILE_RES,
     &             FILE_ROT,FILE_OPT
     
	character*19 FILE_OPTR
	INTEGER FDAT,FDCC,FDEX,FRES,FROT,FOPTR
	
	integer iterp,direcao,rotina
	real resid
	double precision p(npar_max),s(npar_max)

	DIMENSION GF(3*NX/5,3*NX/5),HF(3*NX/5,3*NX/5)
	DIMENSION XF(3*NX/5),YF(3*NX/5),INCF(3*NX/5,2)
	DIMENSION G(NX,NX),H(NX,NX),Hzero(NX,NX)
	DIMENSION X(NX),Y(NX),XM(NX),YM(NX),FI(NX),DFI(NX)
	DIMENSION INC(NX,2),KODE(NX),IT(NX)
	DIMENSION RCONT(NRCX),RCEX(NRCX),fvec(m)
	DIMENSION cx(NPIX),cy(NPIX)
	DIMENSION XIS(NX),U(NX,NBX+1),HBAR(NBX+1,NBX+1),EBAR(NBX+1)
	DIMENSION EEBAR(NBX+1),YIP(NBX+1),CE(NBX+1),ES(NBX+1),XISA(NX)

	DIMENSION DB(NX,NC2MX),NP(NX)
	INTEGER IDCASO(NCasosX,3),I_ELET(16,NELEX)

	DIMENSION NNINC(ninc_max),CINC(ninc_max)
	integer NUMPI(ninc_max), NUMPCI(ninc_max)
	integer VARIA(ninc_max)
        DIMENSION NLD(ninc_max+1,2),INCD(NELDX,3)!para calc do potencial e fluxo nos pts internos
        dimension potno(nx),dpotno(nx,2)!para calc do potencial e fluxo nos pts internos

	double precision dudQ(NX*Nparg_max)
        dimension DFIaux(NX)
	double precision dudQOT(ldfjac,numpar)

	real*8 PAR(numpar),PARG(npar_max)

	integer PCnnp(nparg_max) ! conta qtos nohs parametros cada ponto de controle influencia
	integer PCnp(nparg_max,NX) ! guarda quais nohs parametros cada ponto de controle influencia
	double precision PCgama(nparg_max,NX) ! guarda os pesos com que cada coord dos PC influencia cada noh parametro

	double precision PCgamaSx(nparg_max,NX),PCgamaSy(nparg_max,NX)! guarda os pesos com que cada parametro S dos PCs influencia cada Coordenada X ou Y dos nohs parametros

	integer jNP,cjNPx,cjNPy
	double precision DB2(NX,NC2MX)
	double precision DFI0(NX),DFI1(NX),FI1(NX)
	DIMENSION SOL(NPIX),DSOL(NPIX,2)
	DIMENSION G1(NX,NX),H1(NX,NX)
	logical viavel
	double precision fvecsempena(NRCX)

        real tarray(2),t1sist,t2sist,t1grad,t2grad


        ! Para a rotina de SVD
	real*8 A(NX,NX),Sigma(NX),WORK(5*NX),epssvd
	real*8 Usvd(NX,NX),VT(NX,NX),Vsvd(NX,NX)
	CHARACTER*1 JOBU,JOBV
        real*8 A2(NX,NX),Ainv(NX,NX)

	! para a verificacao da sol do sistema
	real*8 bsist(NX),vRsist(NX)


        double precision Haux(NX,NX),Eaux(NX)
        integer IPIV(NX)

	COMMON /narq/ FDAT,FDCC,FDEX,FRES,FROT
	COMMON /cFDAT/ FILE_DAT,FILE_DEX,FILE_OPT,FILE_RES
   
	COMMON /PG/ GI(8,2),OME(8,2)
	common/plota/ x,y,inc,NN,nfixaux,i_elet,nepe
	common/idr/iterp,direcao,rotina

	COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA
	common/geom/p,s,numpc
	common/dadossensib/PCnnp,PCnp,PCgama,PCgamaSx,PCgamaSy


	save iaval,GF,HF,XF,YF,XM,YM,RCEX,FI,DFI,INCF,NFIX
	save VCC,NCasos,IDCASO,ITCC,NG
	save fvecsempena,inviacont

	integer count1,count2,count_rate,count_max,r



	DATA
     *   (GI(J,1),J=1,4)/0.86113631,-0.86113631,0.33998104,-0.33998104/,
     *   (OME(J,1),J=1,4)/0.34785485,0.34785485,0.65214515,0.65214515/,
     *   (GI(J,2),J=1,8)/0.96028986,-0.96028986,0.79666648,-0.79666648,
     *                   0.52553210,-0.52553210,0.18343464,-0.18343464/,
     *   (OME(J,2),J=1,8)/0.10122854,0.10122854,0.22238103,0.22238103,
     *                    0.31370665,0.31370665,0.36268378,0.36268378/

	data iaval/1/
	data iplota/0/

	! Constante alterando a solucao fundamental
	ALFA = 100.0 
	
	PI= 4.D0*ATAN(1.D0)
	B=  0.D0
	xl_el = .50

	! Parametros para solucao do sistema.

	ITRI= 0
	! WRITE(FRES,110) 
 110  	FORMAT(//,80('*'),/,80('*'),//)


	! Informacoes sobre o numero da entrada de dados,e o numero da
	! saida de resultados.

	FDAT= 12
	FDCC= 13
	FDEX= 14
	FRES= 16
	FROT= 17
	FOPTR=19


	! Monta nomes do arquivo de saida.    
 
	FILE_B   = FILE_DAT
	FILE_DCC = FILE_B//'.dcc'
	
	FILE_B   = FILE_OPT 
	FILE_RES = FILE_B//'.res'
	
	FILE_B   =FILE_OPT
	FILE_OPTR=FILE_B//'_OPT.res'
	
	
	OPEN (FDAT,FILE=FILE_DAT,STATUS='OLD')
	OPEN (FDCC,FILE=FILE_DCC,STATUS='OLD')
	OPEN (FDEX,FILE=FILE_DEX,STATUS='OLD')
	OPEN (FRES,FILE=FILE_RES,STATUS='UNKNOWN')
	open (FOPTR,FILE=FILE_OPTR,status='unknown')


	! INPUT
	IF(IAVAL.EQ.1) THEN
            ! Le os dados globais para as analises
            FILE_ROT = FILE_B//'.rot'
            OPEN (FROT,FILE=FILE_ROT,STATUS='UNKNOWN')

            CALL INPUT0(X,Y,RCEX,VCC,INC,NFIX,NEPE,I_ELET,
     &              NCasos,IDCASO,ITCC,NG,NE)

 125  	    format(5x,'Parametro da solucao fundamental:',e12.5)

            DO I=1,NFIX
                XF(I)= X(I)
                YF(I)= Y(I)
                INCF(I,1)= INC(I,1)
                INCF(I,2)= INC(I,2)
            ENDDO
	ELSE
            DO I=1,NFIX
                X(I)= XF(I)
                Y(I)= YF(I)
                INC(I,1)= INCF(I,1)
                INC(I,2)= INCF(I,2)
            ENDDO
	ENDIF


	if(numpar .eq. numpc)then
            do i=1,numpc
                parg(i*3-2) = p(i*2-1)
                parg(i*3-1) = p(i*2  )
                parg(i*3  ) = (atan(par(i)))/3.1416+0.5
            enddo
	else if(numpar .eq. 2*numpc)then
            do i=1,numpc
                parg(i*3-2) = par(i*2-1)
                parg(i*3-1) = par(i*2  )
                parg(i*3  ) = s(i)
            enddo
	else if(numpar .eq. 3*numpc)then
            do i=1,numpc
                parg(i*3-2) = par(i*3-2)
                parg(i*3-1) = par(i*3-1)
                parg(i*3  ) = (atan(par(i*3)))/3.1416+0.5
            enddo
	endif



	if(iflag .eq. 1 )then
		viavel = .true.

! 		call viabilidade(PARG,X,Y,numpar,numpc,nfix,viavel,xl_el)

                if(iaval .eq. 1 .and. .not. viavel)then
			write(*,*)'Aproximacao inicial nao viavel'
			write(4,*)'Aproximacao inicial nao viavel'
			stop
		endif


		if(viavel)then
			pena = 1.
		else

			write(4,*)'INVIAVEL pena: ',1.+pena
			write(*,*)'INVIAVEL pena: ',1.+pena
	
			do i=1,m
			fvec(i)=fvecsempena(i)*(1.+pena) !penalizacao
			enddo
			iaval = iaval+1
			pena = pena/2.
			return
	    	endif
	endif

!         write(*,*)'parg: ',(PARG(i),i=1,numpc*3)
	CALL GERACONT(PARG,XL_EL,X,Y,NUMPC,FRES,INC,NFIX,NN,N_eq,
     &					npa0) 
     
! ! ! ! ! !   Gera a Geometria no formato do GMSH, possibilitando a criacao de malha do dominio   
!          call togmsh(x,y,xl,inc,NN,FILE_OPT)
! 	  stop
! ! ! ! ! ! ! ! ! ! ! 	 
	 
	 
	nfixaux = NFIX
	write(frot,*) '      Avaliacao:',iaval


        if(iflag.eq.0)then !só plotagem
 	    write(4,*) '                PLOTAGEM No. ',iplota
            resid = 0.
            do i=1,m
                resid = resid + fvec(i)**2
            enddo
            resid = resid/2.0
                
            write(4,222)(PARG(i),i=1,3*numpc)
            write(4,*) 'Residuo:',RESID
            write(4,*)


            call plota_cont(PARG,resid,iplota,numpar,numpc)
            iplota = iplota +1
            return
        endif

        DO IC= 1,NCasos
            ! Introduz as condicoes de contorno com base no caso de carregamento
            WRITE(FRES,1110) 
 1110       FORMAT(//,80('*'),/,80('*'),//)
            WRITE(FRES,1120) IC
 1120       FORMAT(5X,'CASO DE CARREGAMENTO: ',I3)

            CALL DETCC(FI,DFI,VCC,FRES,NEPE,I_ELET,IC,IDCASO,NFIX,KODE)

            ! Formar o sistema de equacoes.
            CALL FMAT(G,H,GF,HF,X,Y,XM,YM,FI,DFI,B,alfa,
     &                KODE,INC,NG,nfix,NN,N_eq,iaval)


!            !!Solucoes do sistema de equacoes.
! 
!                 ! 1) Eliminacao de gauss
!                 call system_clock(count1,count_rate,count_max)
!                 ITRI=0
!                 nc2m = 1
!                 CALL SLNPDM(H,DFI,etol,det,IT,itri,N_eq,nc2m)
!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'SLNPDM(s):',r,'/rate',' NRHS: ',
!      &                  1,' Eqs: ',N_eq 
! 		raux=0.
! 		do i=1,N_eq
! 			raux=raux+dfi(i)**2
! 		enddo
! 		write(*,*)'norma sol: ',dsqrt(raux)
! 		stop
!                 ! fim do solve com eliminacao de gauss


!                 ! 2) GMRES
!                 call system_clock(count1,count_rate,count_max)
!                 NB=    NBX
!                 LMAX=  3
!                 ICONV= 0
!                 IRP=   0
!                 CALL GMRES(H,DFI,XIS,XISA,U,EBAR,EEBAR,HBAR,YIP,CE,
!      &            ES,ETOL,N_eq,NB,NB,LMAX,ICONV,ITER,IRP,FRES)
!                 IF(ICONV.EQ.0) STOP'CONVERGENCIA NAO ATINGIDA - gmres'
!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'GMRES(s):',r,'/rate',' NRHS: ',
!      &                  1,' Eqs: ',N_eq 
!                 !fim do solve com GMRES



                !3) Solve com subrotina do lapack
!                 call system_clock(count1,count_rate,count_max)

                NRHS = 1
                call DGESV(N_eq,NRHS,H,N_eq,IPIV,DFI,N_eq,INFOsolv)

! 			raux=0.
! 			do i=1,N_eq
! 				raux=raux+dfi(i)**2
! 			enddo
! 			write(*,*)'norma sol: ',dsqrt(raux)
!  
!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'DGESV(s):',r,'/rate',' NRHS: ',
!      &                1,' Eqs: ',N_eq 

               ! fim do solve com lapack



C! ! ! ! ! ! 		! processo iterativo para melhor a sol de sistemas com matrizes mal cond (lu- lapack)(GOLUB)
C! ! ! ! ! !		! lembrar de fazer copia da matriz e do vetor dfi antes de resolver o sistema 
C! ! ! ! ! ! 		nvezes = 10
C! ! ! ! ! ! 		do k=1,nvezes
C! ! ! ! ! ! 			
C! ! ! ! ! ! 			call matvec(A,DFI,vRsist,N_eq)
C! ! ! ! ! ! 			do i=1,N_eq
C! ! ! ! ! ! 				vRsist(i) = bsist(i) - vRsist(i)
C! ! ! ! ! ! 			enddo
C! ! ! ! ! ! 
C! ! ! ! ! ! 			raux=0.
C! ! ! ! ! ! 			do i=1,N_eq
C! ! ! ! ! ! 				raux=raux+vRsist(i)**2
C! ! ! ! ! ! 			enddo
C! ! ! ! ! ! 			write(*,*)'    norma do resid da sol: ',dsqrt(raux)
C! ! ! ! ! ! 
C! ! ! ! ! ! 			CALL DGETRS( 'No transpose', N_eq,NRHS,H,N_eq,IPIV, 
C! ! ! ! ! !      &                			vRsist,N_eq,INFOsolve )
C! ! ! ! ! ! 
C! ! ! ! ! ! 			do i=1,N_Eq
C! ! ! ! ! ! 				dfi(i) = dfi(i) + vRsist(i)
C! ! ! ! ! ! 			enddo 
C! ! ! ! ! ! 
C! ! ! ! ! ! 			raux=0.
C! ! ! ! ! ! 			do i=1,N_eq
C! ! ! ! ! ! 				raux=raux+dfi(i)**2
C! ! ! ! ! ! 			enddo
C! ! ! ! ! ! 			write(*,*)'norma sol: ',dsqrt(raux)
C! ! ! ! ! ! 
C! ! ! ! ! ! 		enddo
C! ! ! ! ! ! 		! fim do processo que melhoraria a solucao(Golub)



!              ! 4) Solve com a pseudoinversa para teste:
!                ! verifica a matriz de coeficientes com svd
! 
!                 do i=1,NX
!                     do j=1,NX
!                         A(i,j)=H(i,j)
!                     enddo
!                enddo
! 
! 		epssvd = 0.000001 !valores de sigma da SVD abaix desta tol sao tratados como zero
! 
! 		CALL MPINV (N_eq,N_eq,N_eq,N_eq,A,AINV,
!      &			Sigma,Eaux,Usvd,Vsvd,WORK,IERR,epssvd)
! 
!                 call solvepseudo(Ainv,dfi,N_eq,1)
! 
! 		write(*,*)'Sigma: ',(Sigma(i),i=1,N_eq)
! 		stop
!                ! fim do solve com pseudoinversa da SVD




!  pos processamento







 
            if (metodo.eq.2 .and. iflag.eq.2)then



!                 call system_clock(count1,count_rate,count_max)


                NNP = NN - NFIX !todos os nohs das splines internas sao nohs parametro
			
                do inp=1,NNP
                        NP(inp) = NFIX + inp ! numeros dos nohs parametros
                enddo

! 		Formar o 2o membro do problema das sensibilidades.
                CALL FMATD(X,Y,XM,YM,DB,FI,DFI,alfa,KODE,INC,NG,NNP,NP,
     &						NFIX,NN,N_eq)


                call calcDuDQ(dudQ,DB,PAR,numpc,numpar,N_eq,NFIX,NNP)


!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'montagem(s):',r,'/rate'





!                 !Solve com a pseudoinversa para teste:
!                 call solvepseudo(Ainv,dudq,N_eq,numpar)



!                 ! 1) Eliminacao de gauss
!                 call system_clock(count1,count_rate,count_max)
! 
!                 ITRI=1 !matriz jah triangularizada
!                 nc2m = numpar ! numero de sistemas a resolver
!                 CALL SLNPDM(H,dudQ,etol,det,IT,itri,N_eq,nc2m)
! 
!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'SLNPDM(s):',r,'/rate',' NRHS: ',
!      &                  numpar,' Eqs: ',N_eq 
! 
!                 ! fim do solve com eliminacao de gauss


!                 !2) GMRES
!                 call system_clock(count1,count_rate,count_max)

!                 call preGMRES(H,dudQ,N_eq,numpar,ITER,FRES)

!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'GMRES(s):',r,'/rate',' NRHS: ',
!      &                      numpar,' Eqs: ',N_eq 
!                 ! fim do solve com GMRES    



                !3) Solve com subrotina do lapack - para matria H já triangularizada
!                 call system_clock(count1,count_rate,count_max)

                NRHS = numpar
                CALL DGETRS( 'No transpose', N_eq,NRHS,H,N_eq,IPIV, 
     &                dudQ,N_eq,INFOsolve )

!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'DGESV(s):',r,'/rate',' NRHS: ',
!      &                  NRHS,' Eqs: ',N_eq 

                !fim do solve com lapack




                call sensib_OUTPT_OT(dudQ,dudQOT,numpar,numpc,
     &               ldfjac,I_ELET,IDCASO,IC,NEPE,NCasos,N_eq)


!                 write(*,*)'tempo p calc. do grad.:',(t2grad-t1grad)

            endif !iflag =2

            ! SAIDA DE DADOS
            CALL OUTPT(XM,YM,FI,DFI,KODE,NFIX,NN,N_eq)	

            ! Saida de dados para a otimizacao.
            CALL OUTPT_OT(FI,DFI,RCONT,RCEX,RESID,I_ELET,
     &                IDCASO,IC,NEPE,NCasos,FDEX,FROT)
     
! ! Entrada/ Saida do pós-processamento para o GMSH     
!  	    call INPUTMSH(CX,CY,INCD,NN,NPI,file_b,NLD,NelD,
!      &			nu_msc,nu_msd,nu_trip,nu_trif,ic)    
! ! 	    write(*,*)'pontos internos:',npi
! ! 	    do i=1,npi
! ! 		write(*,*)i,cx(i),cy(i)
! ! 	    enddo
! 	    
! 	    CALL INTER(FI,DFI,CX,CY,X,Y,CINC,SOL,DSOL,KODE,INC,NNINC,
!      &                 NLD,NG,NPI,NN,Nfix,NINCL,C0)
! 
!             CALL OUTPUTMSH(X,Y,XM,YM,FI,DFI,potno,dpotno,CX,CY,SOL,DSOL,
!      &             INC,NNINC,NLD, INCD,NelD,NPI,C0,Nfix,NINCL,N_Eq,NN,
!      &             L,LEC,file_b,nu_msc,nu_msd,nu_trip,nu_trif)
! 
! ! Fim da Entrada/ Saida do pós-processamento para o GMSH       
     
    
        ENDDO  !fim do loop sobre os casos de solicitacao

! 	WRITE(*,1010) ((RCONT(I)),I=1,NRC)
 1010   FORMAT(13E12.4,/,7(13E12.4,/))

	
!  calculo do vetor residuo
        if(metodo.eq.1  .or. (metodo.eq.2 .and. iflag.eq.1))then
            do i=1,m
                fvec(i)=RCONT(i)-RCEX(i)
                fvecsempena(i) = fvec(i)
            enddo
        endif


!  saida 
	if(iflag.eq.1 .or.iflag .eq.2)then
	    write(4,*) '                        Avaliacao:',iaval
            if(numpar .eq. 2*numpc)then
                write(4,222)(PARG(i),i=1,3*numpc)
            else if(numpar .eq. 3*numpc )then
                do i=1,numpc
                    write(4,224)parg(i*3-2),parg(i*3-1),parg(i*3),
     &                                  par(i*3)
                enddo
            elseif(numpar .eq. numpc)then
                do i=1,numpc
                    write(4,224)parg(i*3-2),parg(i*3-1),parg(i*3),
     &                                  par(i)
                enddo
            endif
            write(4,*) 'Residuo:',RESID

            if(iaval.eq.1) write(foptr,1999)
            iaval = iaval+1
            write(FOPTR,*) iaval,resid,(PARG(i),i=1,3*numpc)
        endif
	
! 	call plota_cont(PARG,resid,iplota,numpar,numpc)
!  	stop

 5    format(54f10.5)
 222  format('                      ',3f10.4)
 224  format('                      ',3f10.4,f20.4)
 1999 format('#   Aval.    Resid.       It.Pow.   Direcao   Rotina   ',
     &       ' Num.Eq.   It.GMRES')
 2000 format(i5,5x,f10.6,5i10)

      CLOSE (FDAT)
      CLOSE (FDCC)
      CLOSE (FDEX)
!       CLOSE (FRES)


      RETURN
      END
